/**
 * @file CPEGaussEliminationMethod.h
 * @brief The None principal component Gauss elimination method. 
 * @author XDDDD
 * @version 
 * @date 2021-4-6
 */
#include <iostream>
#include <cmath>

class RationalNum{
private:
		int __a = 0;
		int __b = 1;
public:
		RationalNum(){
			__a = 0;
			__b = 1;
		};
		
		RationalNum(int _a, int _b);

		int get_a() const;
	
		int get_b() const;
		
		void set(int _a, int _b);

		const RationalNum operator=(const RationalNum x);

		friend const RationalNum abs(RationalNum _x);

		friend std::ostream& operator<<(std::ostream &_os, const RationalNum &_t);

		friend RationalNum operator+(const RationalNum _x1, const RationalNum _x2);

		friend RationalNum operator-(const RationalNum _x1, const RationalNum _x2);

		friend RationalNum operator*(const RationalNum _x1, const RationalNum _x2);

		friend RationalNum operator*(int _a, const RationalNum _x); 

		friend RationalNum operator/(const RationalNum _x1, const RationalNum _x2);

		friend bool operator<(const RationalNum _x1, const RationalNum _x2);
		
		friend bool operator>(const RationalNum _x1, const RationalNum _x2);
};

RationalNum::RationalNum(int _a, int _b) {
		int p, q;
		int t;
		p = _a;
		q = _b;
		if(_b == 0) {
			std::cerr << "_b cannot be zero" << std::endl;
			std::exit(-1);
		}
		else if (_b < 0) {
			p = -_a;
			q = -_b;
		}
		if(_a == 0) {
			__a = 0;
			__b = 1;
		}
		else {
			if(p < 0) {
				p = -p;
			}
			while(p % q != 0) {
				t = p % q;
				p = q;
				q = t;
			}
		}
		if(_b < 0) {
			__a = -_a/q;
			__b = -_b/q;
		}
		else {
			__a = _a/q;
			__b = _b/q;
		}
};

int RationalNum::get_a() const {
	return __a;
};

int RationalNum::get_b() const {
	return __b;
};

void RationalNum::set(int _a, int _b) {
	RationalNum t(_a, _b);
	__a = t.__a;
	__b = t.__b;
};

const RationalNum RationalNum::operator=(const RationalNum _x) {
	if(this != &_x) {
		__a = _x.__a;
		__b = _x.__b;
	}
	return *this;
};

std::ostream& operator<<(std::ostream &_os, const RationalNum &_t) {
	if(_t.__b == 1) {
		_os << _t.__a;
	}
	else {
		_os << _t.__a << '/' << _t.__b;
	}
	return _os;
};

RationalNum operator+(const RationalNum _x1, const RationalNum _x2) {
	RationalNum t(_x1.__a*_x2.__b + _x2.__a*_x1.__b, _x1.__b*_x2.__b);
	return t;
};

RationalNum operator-(const RationalNum _x1, const RationalNum _x2) {
	RationalNum t(_x1.__a*_x2.__b - _x2.__a*_x1.__b, _x1.__b*_x2.__b);
	return t;
};

RationalNum operator*(const RationalNum _x1, const RationalNum _x2) {
	RationalNum t(_x1.__a*_x2.__a, _x1.__b*_x2.__b);
	return t;
};

RationalNum operator*(int _a, const RationalNum _x) {
	RationalNum t(_a*_x.__a, _x.__b);
	return t;
}; 

RationalNum operator/(const RationalNum _x1, const RationalNum _x2) {
	if(_x2.__a == 0) {
		std::cerr << "The divisor cannot be zero" << std::endl;
		std::exit(-1);
	}
	RationalNum t(_x1.__a*_x2.__b, _x1.__b*_x2.__a);
	return t;
};
		
bool operator<(const RationalNum _x1, const RationalNum _x2) {
	RationalNum t = _x1 - _x2;
	if(t.__a < 0) {
		return true;
	}
	else {
		return false;
	}
};

bool operator>(const RationalNum _x1, const RationalNum _x2) {
	RationalNum t = _x1 - _x2;
	if(t.__a > 0) {
		return true;
	}
	else {
		return false;
	}
};

const RationalNum abs(RationalNum _x){
	RationalNum t(abs(_x.__a), _x.__b);
	return t;
};

//The input matrix should be stored by row major. A should be a N*N matrix and B should be N*1 matrix.
int CPEGaussEliminationMethod (RationalNum *A, RationalNum *B, int N) {
	int max;
	RationalNum mid;
	for(int j = 1; j != N + 1; j++) {
		max = j;
		for(int t = j + 1; t != N + 1; t++) {
			if(abs(A[t*N - N + j - 1]) > abs(A[j*N - N + j - 1])) {
				max = t;
			}
		}
		if(max != j) {
			for(int t = j; t != N + 1; t++) {
				mid = A[max*N - N + t - 1];
				A[max*N - N + t - 1] = A[j*N - N + t - 1];
				A[j*N - N + t - 1] = mid;
			}
			mid = B[max - 1];
			B[max - 1] = B[j - 1];
			B[j - 1] = mid;
		}
		if(A[j*N - N + j - 1].get_a() == 0) {
			std::cout << "The matrix A is a singularity. The unique solution doest not exist." << std::endl;
			break;
		}
		for(int i = 1; i != N + 1; i++) {
			if(i == j) {
				continue;
			}
			B[i - 1] = B[i - 1] - B[j - 1]*A[i*N - N + j - 1]/A[j*N - N + j - 1];
			for(int k = N; k != j - 1; k--) {
				A[i*N - N + k - 1] = A[i*N - N + k - 1] - A[j*N - N + k - 1]*A[i*N - N + j - 1]/A[j*N - N + j - 1]; 
			}
		}
		
	}
	for(int j = N; j != 1; j--) {
		B[j - 1] = B[j - 1]/A[j*N - N + j - 1];
	}
	B[0] = B[0]/A[0];
	return 0;
};
